Bloque 1: Modelos de tidymodels

Analítica basada en árboles de clasificación y regresión

Modelo

El objetivo de un modelo es proporcionar un resumen simple de baja dimensión de un conjunto de datos. Idealmente, el modelo capturará patrones generados por el fenómeno de interés e ignorará el “ruido” (es decir, la variación aleatoria que no le interesa).

Tradicionalmente, el enfoque del modelado está en la inferencia o confirmar que una hipótesis es cierta. Hay un par de ideas que nos pueden ayudar a hacer inferencias correctamente:

1. Cada observación puede usarse para exploración o confirmación, no ambas.

2. Puede utilizar una observación tantas veces como desee para la exploración, pero solo puede utilizarla una vez para la confirmación. Tan pronto como utilice una observación dos veces, pasará de la confirmación a la exploración.

Esto es necesario porque para confirmar una hipótesis debe utilizar datos independientes de los datos que utilizó para generar la hipótesis. De lo contrario, será demasiado optimista. No hay absolutamente nada de malo en la exploración, pero nunca debe vender un análisis exploratorio como un análisis confirmatorio porque es fundamentalmente engañoso.

Si realmente desea realizar un análisis confirmatorio, un enfoque es dividir sus datos en tres partes antes de comenzar el análisis:

1. El 60% de sus datos se destina a un conjunto de entrenamiento, training, (o exploración). Puede hacer lo que quiera con estos datos: visualizarlos y ajustarles toneladas de modelos.

2. El 20% se destina a un conjunto de consultas, query. Puede usar estos datos para comparar modelos o visualizaciones a mano, pero no puede usarlos como parte de un proceso automatizado.

3. El 20% se retiene para una prueba testing. Solo puede usar estos datos UNA VEZ, para probar su modelo final.

Esta partición le permite explorar los datos de entrenamiento, generando ocasionalmente hipótesis candidatas que verifica con el conjunto de consultas. Cuando esté seguro de que tiene el modelo correcto, puede verificarlo una vez con los datos de prueba.

Tidymodels

Tidymodels, es una interfaz que unifica bajo un único marco cientos de funciones de distintos paquetes, facilitando en gran medida todas las etapas de preprocesado, entrenamiento, optimización y validación de modelos predictivos. Los paquetes principales que forman parte del ecosistema tidymodels son:

• parsnip: para la definición de modelos. Implementación tidy del caret.

• recipes: para el preprocesado de datos y feature engineering.

• rsample: para validar los modelos por métodos de resampling.

• dials: para crear y manejar el valor de los hiperparámetros.

• tune: para hacer tuning de modelos.

• yardstick: para calcular métricas de modelos.

• workflows: para combinar todos los pasos del preprocesado y modelado en un único objeto.

El procesamiento lo realizan los paquetes rsample y recipes. El modelo principal está en parsnip (el equivalente de Tidy a caret) y para la validación Yardstick. Cabe señalar que caret sigue teniendo la mejor implementación de las tablas de matrices de confusión, con un paquete adaptado ConfusionTableR.

División Train y test

Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo se acercan las predicciones a los valores verdaderos de la variable respuesta. Para poder cuantificarlo de forma correcta, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya “visto”, es decir, que no hayan participado en su ajuste. Con esta finalidad, se divide el conjunto de datos disponibles en un conjunto de entrenamiento (train) y un conjunto de prueba (test).

El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80%-20% suele dar buenos resultados. El reparto debe hacerse de forma aleatoria o aleatoria-estratificada. Evaluación de modelos (ver link).

Para este ejemplo, analizaremos el conjunto de datos SaratogaHiyses, que se encuentra en el paquete nisaucData, el cual contiene información sobre los precio de 1728 viviendas de Saratoga County, New York, USA en el año 2006. Además del precio, este dataset incluye otras 15 variables:

• price: precio de la vivienda.

• lotSize: metros cuadrados de la vivienda.

• age: antigüedad de la vivienda.

• landValue: valor del terreno.

• livingArea: metros cuadrados habitables.

• pctCollege: porcentaje del vecindario con título universitario.

• bedrooms: número de dormitorios.

• firplaces: número de chimeneas.

• bathrooms: número de cuartos de baño (el valor 0.5 hace referencia a cuartos de baño sin ducha).

• rooms: número de habitaciones.

• heating: tipo de calefacción.

• fuel: tipo de alimentación de la calefacción (gas, electricidad o diesel).

• sewer: tipo de desagüe.

• waterfront: si la vivienda tiene vistas al lago.

• newConstruction: si la vivienda es de nueva construcción.

• centralAir: si la vivienda tiene aire acondicionado.

El objetivo es obtener un modelo capaz de predecir el precio del alquiler casa. Por lo tanto para poder llamar los datos tomamos el siguiente código:

data("SaratogaHouses", package = "mosaicData")
data <- SaratogaHouses
colnames(data)

 [1] "price"           "lotSize"         "age"             "landValue"       "livingArea"     
 [6] "pctCollege"      "bedrooms"        "fireplaces"      "bathrooms"       "rooms"          
[11] "heating"         "fuel"            "sewer"           "waterfront"      "newConstruction"
[16] "centralAir"     

Con este conjunto de datos procedemos a hacer la división entre train y test.

library(tidymodels)

Registered S3 method overwritten by 'tune':
  method                   from   
  required_pkgs.model_spec parsnip
── Attaching packages ────────────────────────────────────────────────────────────── tidymodels 0.1.3 ──
✓ broom        0.7.8      ✓ recipes      0.1.16
✓ dials        0.0.9      ✓ rsample      0.1.0 
✓ dplyr        1.0.7      ✓ tibble       3.1.2 
✓ ggplot2      3.3.5      ✓ tidyr        1.1.3 
✓ infer        0.5.4      ✓ tune         0.1.5 
✓ modeldata    0.1.0      ✓ workflows    0.2.2 
✓ parsnip      0.1.6      ✓ workflowsets 0.0.2 
✓ purrr        0.3.4      ✓ yardstick    0.0.8 
── Conflicts ───────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
x purrr::discard() masks scales::discard()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x recipes::step()  masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.

# Reparto de datos en train y test
set.seed(123)
split_inicial <- initial_split(
                    data   = data, # datos que queremos dividir
                    prop   = 0.8, # Proporción del 80%
                    strata = price # Variable respuesta 
                 )
data_train <- training(split_inicial)
data_test  <- testing(split_inicial)

Es importante verificar que la distribución de la variable respuesta es similar en el conjunto de entrenamiento y en el de test. Para asegurar que esto se cumple, la función initial_split() permite identificar con el argumento strata la variable en base a la cual hacer el reparto.

summary(data_train$price)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5000  145000  189050  212423  259000  775000 

summary(data_test$price)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  20000  145000  189950  210157  257750  620000 

Este tipo de reparto estratificado asegura que el conjunto de entrenamiento y el de test sean similares en cuanto a la variable respuesta, sin embargo, no garantiza que ocurra lo mismo con los predictores. Por ejemplo, en un set de datos con 100 observaciones, un predictor binario que tenga 90 observaciones de un grupo y solo 10 de otro, tiene un alto riesgo de que, en alguna de las particiones, el grupo minoritario no tenga representantes. Si esto ocurre en el conjunto de entrenamiento, algunos algoritmos darán error al aplicarlos al conjunto de test, ya que no entenderán el valor que se les está pasando. Este problema puede evitarse eliminando variables con varianza próxima a cero, lo cual veremos más adelante.

Exclusión de variables con varianza próxima a cero

No se deben incluir en el modelo predictores que contengan un único valor (varianza-cero) ya que no aportan información. Tampoco es conveniente incluir predictores que tengan una varianza próxima a cero, es decir, predictores que toman solo unos pocos valores, de los cuales, algunos aparecen con muy poca frecuencia. El problema con estos últimos es que pueden convertirse en predictores con varianza cero cuando se dividen las observaciones por validación cruzada o bootstrap.

¿Recuerdas que lo vimos en el análisis exploratorio, EDA?

Aquí, también podemos utilizar la función step_nzv() del paquete recipes que identifica predictores potencialmente problemáticos, es decir, aquellos que tiene un único valor o varianza- cero que cumplen dos condiciones:

• Variación de frecuencia: Es la variación entre la frecuencia del valor más común y la frecuencia del según valor más común. Este ratio tiende a 1 si las frecuencias están equidistribuidas y a valores grandes cuando la refuencia del valor mayoritario supera por mucho al resto (el denominador es un número decimal pequeño). Valor por defecto freq_cut = 95/5.

• Porcentaje de valores únicos: número de valores únicos divido entre el total de la muestra (multiplicado por 100). Esta valor se acerca a cero cuando la varianza-cero es mayor. Valor por defecto uniqueCut = 10.

Si bien, cuando eliminamos predictores no informativos podría considerarse un paso propio del proceso de selección de predictores, dado que consiste en un filtrado por varianza, tiene que realizarse antes de estandarizar los datos, ya que después, todos los predictores tienen varianza 1.

Estandarización y escalado de variables numéricas

Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Muchos algoritmos de machine learning (SVM, redes neuronales, lasso, etc) son sensibles a esto, de forma que, si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:

• Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.

• Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala. Hay dos formas de lograrlo:

  • Normalización Z-score: dividir cada predictor entre su desviación típica después de haber sido centrado, de esta forma, los datos pasan a tener una distribución normal.

  \[ z = \frac{x - \mu}{\sigma} \]

  • Estandarización max-min: transformar los datos de forma que estén dentro del rango [0, 1].

\[ x_{nom}=\frac{x-x_{min}}{x_{max}-x_{min}} \] NOTA: Nunca se debe estandarizar las variables después de ser binarizadas.

Binarización de las variables cualitativas

La binarización consiste en crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. A este proceso también se le conoce como one hot encoding. Por ejemplo, una variable llamada color que contenga los niveles rojo, verde y azul, se convertirá en tres nuevas variables (color_rojo, color_verde, color_azul), todas con el valor 0 excepto la que coincide con la observación, que toma el valor 1.

Por defecto, la función step_dummy(all_nominal(), -all_outcomes()) binariza todas las variables almacenadas como tipo factor o character, excepto la variable respuesta. Además, elimina uno de los niveles para evitar redundancias. Volviendo al ejemplo anterior, no es necesario almacenar las tres variables, ya que, si color_rojo y color_verde toman el valor 0, la variable color_azul toma necesariamente el valor 1. Si color_rojo o color_verde toman el valor 1, entonces color_azul es necesariamente 0.

El paquete recipes incorpora una amplia variedad de funciones para preprocesar los datos, facilitando el aprendizaje de las transformaciones únicamente con observaciones de entrenamiento, y poder aplicarlas después a cualquier conjunto de datos. La idea detrás de este paquete es la siguiente:

1. Definir cuál es la variable respuesta, los predictores y el conjunto de datos de entrenamiento, recipe().

2. Definir todas las transformaciones (escalado, selección, filtrado…) que se desea aplicar, step_().

3. Aprender los parámetros necesarios para dichas transformaciones con las observaciones de entrenamiento rep().

4. Aplicar las transformaciones aprendidas a cualquier conjunto de datos juice(), bake().

# Se almacenan en un objeto `recipe` todos los pasos de preprocesado y, finalmente, se aplican a los datos.
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

transformer

Data Recipe

Inputs:

Operations:

Removing rows with NA values in all_predictors()
Sparse, unbalanced variable filter on all_predictors()
Centering for all_numeric(), -all_outcomes()
Scaling for all_numeric(), -all_outcomes()
Dummy variables from all_nominal(), -all_outcomes()

Una vez que se ha definido el objeto recipe, con la función prep() se aprenden las transformaciones con los datos de entrenamiento y se aplican a los dos conjuntos con bake().

# Se entrena el objeto recipe
transformer_fit <- prep(transformer)

# Se aplican las transformaciones al conjunto de entrenamiento y de test
data_train_prep <- bake(transformer_fit, new_data = data_train)
data_test_prep  <- bake(transformer_fit, new_data = data_test)

glimpse(data_train_prep)

Rows: 1,380
Columns: 18
$ lotSize                 <dbl> -0.43894231, 0.26607619, 0.48189818, 0.61139137, 12.19383812, -0.51088…
$ age                     <dbl> 3.62268073, 0.11471419, 0.28667334, -0.91704067, -0.50433872, -0.64190…
$ landValue               <dbl> -0.7674967, -0.5829431, -0.3570715, -0.3543170, -0.8363600, -0.9603140…
$ livingArea              <dbl> 0.31203819, -0.96981699, -0.19293506, -0.20588309, -1.69490679, -0.730…
$ pctCollege              <dbl> -0.4558859, -3.2928560, -0.4558859, -0.4558859, -1.4341515, -1.4341515…
$ bedrooms                <dbl> 1.0465989, 1.0465989, -0.1681248, -0.1681248, -1.3828485, -0.1681248, …
$ fireplaces              <dbl> 0.7295147, 0.7295147, -1.0876880, -1.0876880, -1.0876880, -1.0876880, …
$ bathrooms               <dbl> -1.3602353, -1.3602353, -0.6045490, 0.1511373, -1.3602353, -0.6045490,…
$ rooms                   <dbl> 0.425229547, 0.425229547, 0.425229547, -0.440920675, -1.307070897, 0.4…
$ price                   <int> 109000, 120000, 90000, 120000, 85860, 127000, 89900, 60000, 87500, 112…
$ heating_hot.water.steam <dbl> 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1…
$ heating_electric        <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ fuel_electric           <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ fuel_oil                <dbl> 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ sewer_public.commercial <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1…
$ sewer_none              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ newConstruction_No      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ centralAir_No           <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1…

Tras el preprocesado de los datos, se han generado un total de 18 variables (17 predictores y la variable respuesta).

Modelado

Ya analizado los datos (entenderlos por medio del EDA), que han sido preprocesados y los predictores seleccionados, el siguiente paso es emplear un algoritmo de machine learning que permita crear un modelo capaz de representar los patrones presentes en los datos de entrenamiento y generalizarlos a nuevas observaciones. Encontrar el mejor modelo no es fácil, existen multitud de algoritmos, cada uno con unas características propias y con distintos parámetros que deben ser ajustados. Por lo general, las etapas seguidas para obtener un buen modelo son:

• Ajuste/entrenamiento: consiste en aplicar un algoritmo de machine learning a los datos de entrenamiento para que el modelo aprenda.

• Evaluación/validación: el objetivo de un modelo predictivo no es ser capaz de predecir observaciones que ya se conocen, sino nuevas observaciones que el modelo no ha visto. Para poder estimar el error que comete un modelo es necesario recurrir a estrategias de validación entre las que destacan: validación simple, bootstrap y validación cruzada.

• Optimización de hiperparámetros: muchos algoritmos de machine learning contienen en sus ecuaciones uno o varios parámetros que no se aprenden con los datos, a estos se les conoce como hiperparámetros. Por ejemplo, los SVM lineal tiene el hiperparámetro de coste C y los árboles de regresión la profundidad del árbol tree_depth y el número mínimo de observaciones por nodo min_n. No existe forma de conocer de antemano cuál es el valor exacto de un hiperparámetro que da lugar al mejor modelo, por lo que se tiene que recurrir a estrategias de validación para comparar distintos valores.

• Predicción: una vez creado el modelo, este se emplea para predecir nuevas observaciones.

Es a lo largo de todo este proceso donde más destacan las funcionalidades ofrecidas por tidymodels, permitiendo emplear la misma sintaxis para ajustar, optimizar, evaluar y predecir un amplio abanico de modelos, variando únicamente el nombre del algoritmo.

Aunque tidymodels permite todo esto con apenas unas pocas líneas de código, son muchos los argumentos que pueden ser adaptados, cada uno con múltiples posibilidades. Con el objetivo de exponer mejor cada una de las opciones, en lugar de crear directamente un modelo final, se muestran ejemplos de menor a mayor complejidad.

Entrenamiento

Los modelos de tidymodels son accesibles por medio del paquete parsnip. Aunque un mismo modelo puede estar en varios paquetes con nombres distintos.

Para crear un modelo se requieren únicamente tres pasos:

• Definir el tipo de modelo y sus parámetros

• Definir qué implementación del modelo se quiere aplicar (engine)

• Ajustar el modelo

Vamos a ajustar un modelo de árbol de regresión para predecir el precio de la vivienda de Saratoga en función de todos los predictores disponibles.
A excepción de la implementación del algoritmo ("rpart"), rpart (viene por default) o "C5.0" (solo para clasificación), todos los demás argumentos de la función decision_tree() se dejan por defecto.

model_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")
model_tree

Decision Tree Model Specification (regression)

Computational engine: rpart 

Con las líneas anteriores sólo guardamos el tipo de algoritmo (sus parámetros, hiperparámetros y el paquete en el que está implementado). El sigiente paso es ajustar el modelo, para los cual tenemos dos opciones:

• fit() que emplea como argumento una formula para definir la variable respuesta y los predictores.

• fit_xy() que emplea los argumentos x y y para definir la matriz de predictores y el vector con la variable respuesta.

# Entrenamiento empleando fórmula
model_tree_fit <- model_tree %>%
                   fit(
                     formula = price ~ .,
                     data    = data_train_prep
                   )

# Entrenamiento empleando x e Y.
variable_respuesta <- "price"

predicores <- setdiff(colnames(data_train_prep), variable_respuesta)

model_tree_fit <- model_tree %>%
                   fit_xy(
                     x = data_train_prep[, predicores],
                     y = data_train_prep[[variable_respuesta]]
                   )

Para ver nuestro modelo utilizamos lo que se almacena dentro de model_tree_fit como fit:

model_tree_fit$fit

n= 1380 

node), split, n, deviance, yval
      * denotes terminal node

  1) root 1380 1.382982e+13 212423.0  
    2) livingArea< 0.3856801 956 3.951824e+12 172884.3  
      4) landValue< -0.3405443 540 1.320059e+12 147526.2  
        8) livingArea< -0.3968666 397 5.752973e+11 135085.3 *
        9) livingArea>=-0.3968666 143 5.127270e+11 182065.0 *
      5) landValue>=-0.3405443 416 1.833788e+12 205801.0  
       10) landValue< 0.8136044 359 1.118047e+12 196236.6  
         20) livingArea< -0.001951587 251 7.002854e+11 182517.3 *
         21) livingArea>=-0.001951587 108 2.607222e+11 228121.3 *
       11) landValue>=0.8136044 57 4.760637e+11 266039.8 *
    3) livingArea>=0.3856801 424 5.013749e+12 301571.5  
      6) landValue< 0.9733672 311 1.944443e+12 268743.7  
       12) livingArea< 1.830195 277 1.304977e+12 256947.8  
         24) centralAir_No>=0.5 140 4.264425e+11 233996.8 *
         25) centralAir_No< 0.5 137 7.294296e+11 280401.5 *
       13) livingArea>=1.830195 34 2.869182e+11 364844.9 *
      7) landValue>=0.9733672 113 1.811733e+12 391920.8  
       14) rooms< 1.940992 89 1.063606e+12 365619.8  
         28) age< -0.8998448 41 9.079364e+10 317584.5 *
         29) age>=-0.8998448 48 7.974026e+11 406650.0  
           58) livingArea< 0.7029069 13 7.632773e+10 320607.7 *
           59) livingArea>=0.7029069 35 5.890850e+11 438608.6  
            118) age>=-0.7622774 27 2.569939e+11 403566.7 *
            119) age< -0.7622774 8 1.870413e+11 556875.0 *
       15) rooms>=1.940992 24 4.582587e+11 489453.6 *

Validación del modelo

Lo que buscamos ver es que tan certera es la predicción de nuestro modelo, ya sea con observaciones que no ha visto el modelo y por ende, con observaciones futuras.

El error mostrado por defecto tras entrenar un modelo suele ser el error de entrenamiento, el error que comete el modelo al predecir las observaciones que ya ha “visto”.

Estos errores son útiles para entender cómo está aprendiendo el modelo (estudio de residuos), no es una estimación realista de cómo se comporta el modelo ante nuevas observaciones (el error de entrenamiento suele ser demasiado optimista).

Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. El paquete rsampler incorpora los métodos Bootstrap (bootstraps), V-Fold Cross-Validation y Repeated V-Fold Cross-Validation (vfold_cv), Nested or Double Resampling (nested_cv), Group V-Fold Cross-Validation (group_vfold_cv), Leave-One-Out Cross-Validation (loo_cv), Monte Carlo Cross-Validation(mc_cv). Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea: ajustar y evaluar el modelo de forma repetida, empleando cada vez distintos subconjuntos creados a partir de los datos de entrenamiento y obteniendo en cada repetición una estimación del error (simulaciones de los datos). El promedio de todas las estimaciones tiende a converger en el valor real del error de test.

Se ajusta de nuevo el modelo, esta vez con validación cruzada repetida para estimar su error.

En primer lugar, se crea un objeto resample que contenga la información sobre las observaciones que forman parte de cada partición. Dado que el reparto es aleatorio, es importante emplear una semilla set.seed() para que los resultados sean reproducibles.

set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              repeats = 10,
              strata  = price
            )
head(cv_folds)

Con la función fit_resamples() el modelo se ajusta y evalúa con cada una de las particiones de forma automática, calculando y almacenando en cada iteración la métrica de interés. Como se comentó anteriormente, cuando se realizan validaciones por resampling, el preprocesado debe ocurrir dentro de cada iteración.

Para seguir este principio, las particiones se crean con el conjunto de datos de entrenamiento sin preprocesar y se le pasa al argumento preprocessor un objeto recipe que contenga la definición del preprocesado.

model_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")

validation_fit <- fit_resamples(
                    object       = model_tree,
                    # El objeto recipe no tiene que estar entrenado
                    preprocessor = transformer,
                    # Las resamples se tienen que haber creado con los datos sin 
                    # prerocesar
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse, mae),
                    control      = control_resamples(save_pred = TRUE)
                  )

head(validation_fit)

Los resultados de fit_resamples() se almacenan en forma de tibble, donde las columnas contienen la información sobre cada partición: su id, las observaciones que forman parte, las métricas calculadas, si ha habido algún error o warning durante el ajuste, y las predicciones de validación si se ha indicado control = control_resamples(save_pred = TRUE). La información puede extraerse haciendo un unnest() de las columnas o empleando las funciones auxiliares collect_predictions() y collect_metrics().

# Métricas promedio de todas las particiones
validation_fit %>% collect_metrics(summarize = TRUE)

# Métricas individuales de cada una de las particiones
validation_fit %>% collect_metrics(summarize = FALSE) %>% 
  head()

library("ggpubr")

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

# Valores de validación (mae y rmse) obtenidos en cada partición y repetición.
p1 <- ggplot(
        data = validation_fit %>% collect_metrics(summarize = FALSE),
        aes(x = .estimate, fill = .metric)) +
      geom_density(alpha = 0.5) +
      theme_bw() 
p2 <- ggplot(
        data = validation_fit %>% collect_metrics(summarize = FALSE),
        aes(x = .metric, y = .estimate, fill = .metric, color = .metric)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.1) +
      geom_jitter(width = 0.05, alpha = 0.3) +
      coord_flip() +
      theme_bw() +
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  
ggarrange(p1, p2, nrow = 2, common.legend = TRUE, align = "v") %>% 
annotate_figure(
  top = text_grob("Distribución errores de validación cruzada", size = 15)
)

El rmse (raíz del error medio cuadrático) promedio estimado mediante validación cruzada repetida es de 67990. Este valor será contrastado más adelante cuando se calcule el rmse del modelo con el conjunto de test.

Almacenar las predicciones de cada partición es útil para poder evaluar los residuos del modelo y diagnosticar así su comportamiento.

# Predicciones individuales de cada observación.
# Si summarize = TRUE se agregan todos los valores predichos a nivel de
# observación.
validation_fit %>% collect_predictions(summarize = TRUE) %>% head()

p1 <- ggplot(
        data = validation_fit %>% collect_predictions(summarize = TRUE),
        aes(x = price, y = .pred)
      ) +
      geom_point(alpha = 0.3) +
      geom_abline(slope = 1, intercept = 0, color = "firebrick") +
      labs(title = "Valor predicho vs valor real") +
      theme_bw()


p2 <- ggplot(
        data = validation_fit %>% collect_predictions(summarize = TRUE),
        aes(x = .row, y = price - .pred)
      ) +
      geom_point(alpha = 0.3) +
      geom_hline(yintercept =  0, color = "firebrick") +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(
        data = validation_fit %>% collect_predictions(summarize = TRUE),
        aes(x = price - .pred)
      ) +
      geom_density() + 
      labs(title = "Distribución residuos del modelo") +
      theme_bw()

p4 <- ggplot(
        data = validation_fit %>% collect_predictions(summarize = TRUE),
        aes(sample = price - .pred)
      ) +
      geom_qq() +
     geom_qq_line(color = "firebrick") +
      labs(title = "Q-Q residuos del modelo") +
      theme_bw()

ggarrange(plotlist = list(p1, p2, p3, p4)) %>%
annotate_figure(
  top = text_grob("Distribución residuos", size = 15, face = "bold")
)

Cuando se aplican métodos de* resampling debemos tener en cuenta dos cosas:

1. El coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto,

2. y la reproducibilidad en la creación de las particiones.

fit_resamples() permite paralelizar el proceso si se registra un backed paralelo, para ello usamos la función resterDoParallel() del paquete doParallel.

library("doParallel")

Loading required package: foreach
Warning: package ‘foreach’ was built under R version 4.0.5
Attaching package: ‘foreach’

The following objects are masked from ‘package:purrr’:

    accumulate, when

Loading required package: iterators
Loading required package: parallel

registerDoParallel(cores = detectCores() - 1)
set.seed(2020)
model_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")

validacion_fit <- fit_resamples(
                    object       = model_tree,
                    preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse, mae),
                    control      = control_resamples(save_pred = TRUE)
                  )

stopImplicitCluster()

Hiperparámetros (tuning)

Muchos modelos, entre ellos los árboles de regresión, contienen parámetros (datos de entrada al modelo) que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. ¡Ojo! no se trata de tener un modelo y correrlo y ya con esto llamarnos científicos de datos.

Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades por medio de algoritmos de optimización.

Podemos seguir los siguientes pasos:

1. Escoger un conjunto de valores para el o los hiperparámetros.

2. Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación.

3. Ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.

El modelo decision_tree empleado hasta ahora tienen tres hiperparámetros: cost_complexity, tree_depth y min_n. Todos ellos controlan el tamaño del árbol final. Por ejemplo, cuanto mayor el valor de tree_depth, más ramificaciones tiene el árbol y más flexible es el modelo. Esto le permite ajustarse mejor a las observaciones de entrenamiento pero con el riesgo de overfitting. Con las funciones tune() y tune_grid() del paquete tune se pueden explorar diferentes valores de hiperparámetros sin apenas tener que cambiar la estructura del código.

Se vuelve a ajustar un modelo decision_tree con diferentes valores de tree_depth y min_n, y empleando validación cruzada repetida para identificar con cuáles de ellos se obtienen mejores resultados.

# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = model_tree,
              # El objeto recipe no tiene que estar entrenado
              preprocessor = transformer,
              # Las resamples se tienen que haber creado con los datos sin 
              # prerocesar
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_grid(save_pred = TRUE),
              # Número de combinaciones generadas automáticamente
              grid         = 70
            )
stopImplicitCluster()

Los resultados de la búsqueda de hiperparámetros pueden verse haciendo un unnest() de el tibble que se crea de la ejecución grid_fit, o con las funciones auxiliares collect_metrics(), collect_predictions, show_best() y select_best().

grid_fit %>% unnest(.metrics) %>% head()

grid_fit %>% collect_metrics(summarize = TRUE) %>% head()

grid_fit %>% show_best(metric = "rmse", n = 5)

Es conveniente visualizar la evolución del error en función de los hiperparámetros. Sin embargo, hay que tener en cuenta que estos no son independientes los unos de los otros, el comportamiento de un hiperparámetro puede cambiar mucho dependiendo del valor que tomen los otros.

grid_fit %>%
  collect_metrics(summarize = TRUE) %>%
  filter(.metric == "rmse") %>%
  select(-c(.estimator, n)) %>%
  pivot_longer(
    cols = c(tree_depth, min_n),
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(x = value, y = mean, color = parameter)) +
  geom_point() +
  geom_line() + 
  labs(title = "Evolución del error en función de los hiperparámetros") +
  facet_wrap(facets = vars(parameter), nrow = 2, scales = "free") +
  theme_bw() + 
  theme(legend.position = "none")

grid_fit %>%
  collect_metrics(summarize = TRUE) %>%
  filter(.metric == "rmse") %>%
  select(-c(.estimator, n)) %>%
  ggplot(aes(x = tree_depth, y = min_n, color = mean, size = mean)) +
  geom_point() +
  scale_color_viridis_c() +
  labs(title = "Evolución del error en función de los hiperparámetros") +
  theme_bw()

En este caso, parece que el error del modelo se reduce a medida que el valor de min_n aumenta (aunque no es estable) y también el de tree_depth. En este último la mejora parece estabilizarse a partir de 4.

Aunque los autores de tidymodels han establecido por defecto valores adecuados basados en su amplia experiencia, en muchos casos es conveniente tener más control sobre la búsqueda, ya que la selección se hace de forma automática. Para conseguir este control, podemos utilizar expand.grid() (para especificar exactamente qué valores se emplean), o con las funciónes regular_grid(), grid_random(), grid_max_entropy(), y grid_latin_hypercube() del paquete dials.

Se repite la búsqueda de mejores hiperparámetros pero esta vez definiendo el espacio de búsqueda.

# Definición del modelo y de los hiperparámetros a optimizar
model_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Grid de Hiperparámetros
set.seed(1234)
hiperpar_grid <- grid_random(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range      = c(2, 100), trans = NULL),
                  # Número combinaciones aleatorias probadas
                  size = 50
                )

# Ejecución y optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = model_tree,
              # El objeto recipe no tiene que estar entrenado
              preprocessor = transformer,
              # Las resamples se tienen que haber creado con los datos sin 
              # prerocesar
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid         = hiperpar_grid
            )
stopImplicitCluster()

Ahora vemos los resultados:

grid_fit %>% show_best(metric = "rmse", n = 5)

grid_fit %>%
  collect_metrics(summarize = TRUE) %>%
  filter(.metric == "rmse") %>%
  select(-c(.estimator, n)) %>%
  pivot_longer(
    cols = c(tree_depth, min_n),
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(x = value, y = mean, color = parameter)) +
  geom_point() +
  geom_line() + 
  labs(title = "Evolución del error en función de los hiperparámetros") +
  facet_wrap(facets = vars(parameter), nrow = 2, scales = "free") +
  theme_bw() + 
  theme(legend.position = "none")

Los mejores resultados, en base al rmse se han obtenido con los hiperparámetros:

# Selección de los mejores hiperparámetros encontrados
best_hip <- select_best(grid_fit, metric = "rmse")
best_hip

Modelo final

Ya haciendo la optimización de los hiperparámetros, entramos al modelo pero con los datos de entrenamiento con estos parámetros.

final_model_tree <- finalize_model(x = model_tree, 
                                   parameters = best_hip)
final_model_tree

Decision Tree Model Specification (regression)

Main Arguments:
  tree_depth = 10
  min_n = 69

Computational engine: rpart 

final_model_tree_fit  <- final_model_tree %>%
                          fit(
                            formula = price ~ .,
                            data    = data_train_prep
                          )

Predicción

Una vez que el modelo final ha sido ajustado, con la función predict() se pueden predecir nuevos datos. Los argumentos de la función son:

• object: un modelo entrenado.

• newdata: un dataframe con nuevas observaciones.

• type: tipo de predicción ("numeric", "class", "prob", "conf_int", "pred_int", "quantile", or "raw").

predicciones <- final_model_tree_fit  %>%
                predict(
                  new_data = data_test_prep,
                  #new_data = bake(transformer_fit, datos_test),
                  type = "numeric"
                )
predicciones %>% head()

Error de test

Ya hemos validado con modelos de validación (CV, bootstrapping…) y se han conseguido estimaciones del error del modelo al predecir nuevas observaciones. Pero la mejor forma de evaluar el modelo final es prediciendo un conjunto test, es decir, utilizando un conjunto de datos que no hemos usado en el entrenamiento y en la optimización. El paquete yardstick tiene funciones predefinidas para métricas en cuestión que nos pueden ayudar.

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

rmse(predicciones, truth = price, estimate = .pred, na_rm = TRUE)

En el apartado de optimización, se estimó, mediante validación cruzada repetida, que el rmse del modelo model_tree era de 67790, un valor próximo al obtenido con el conjunto de test baja a 61536.

Workflows

Los workflows permiten combinar en un sólo objeto todos los elementos que se encargan del procesamiento (recipes), modelado (parsnip y tune) y post-procesado. Para crear el workflows se van encadenando los elementos con las funciones add_* o modificando los elementos de las funciones update_*.

# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# Definición de procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(model_tree)

workflow_modelado

══ Workflow ════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_naomit()
• step_nzv()
• step_center()
• step_scale()
• step_dummy()

── Model ───────────────────────────────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (regression)

Main Arguments:
  tree_depth = tune()
  min_n = tune()

Computational engine: rpart 

El objeto workflow puede ajustarse directamente con la función fit() o hacer tuning de hiperparámetros con tune_grid().

# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range = c(2, 100), trans = NULL),
                  # Número valores por hiperparámetro
                  levels = c(3, 3)
                )

# Ejecución optimización hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = workflow_modelado,
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid         = hiperpar_grid
            )
stopImplicitCluster()

# Entrenamiento final
best_hiper <- select_best(grid_fit, metric = "rmse")

final_model_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = best_hiper
                    ) %>%
                    fit(
                      data = data_train
                    ) %>%
                    pull_workflow_fit()

Estrategias de tuning

Uno de los problemas que tenemos para encontrar los mejores hiperparámetros es el coste computacional. Una alternativa es la búsqueda de hiperparámetros con métodos de optimización bayesiana.

En términos generales, la optimización bayesiana de hiperparámetros consiste en crear un modelo probabilístico en el que la función objetivo es la métrica de validación del modelo (rmse, auc, precisión, etc). Con esta estrategia, se consigue que la búsqueda se vaya redirigiendo en cada iteración hacia las regiones de mayor interés.

El objetivo final es reducir el número de combinaciones de hiperparámetros con las que se evalúa el modelo, eligiendo únicamente los mejores candidatos. Esto significa que, la ventaja frente a las otras estrategias mencionadas, se maximiza cuando el espacio de búsqueda es muy amplio o la evaluación del modelo es muy lenta.

# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
                ) %>%
               set_engine(engine = "rpart")

# Procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(model_tree)

# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range = c(2, 100), trans = NULL),
                  # Número valores por hiperparámetro
                  levels = c(3, 3)
                )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_bayes(
              workflow_modelado, 
              resamples = cv_folds,
              # Iniciación aleatoria con 20 candidatos
              initial = 20,
              # Numero de iteraciones de optimización
              iter    = 30,
              # Métrica optimizada
              metrics = metric_set(rmse),
              control = control_bayes(no_improve = 20, verbose = FALSE)
            )

! No improvement for 20 iterations; returning current results.

stopImplicitCluster()

Al igual que con el método anterior, podemos ver la evolución del error:

autoplot(grid_fit, type = "performance") +
  labs(title = "Evolución del error") +
  theme_bw()

Para mostrar los mejores hiperparámetros encontrados

show_best(x = grid_fit, metric = "rmse")

Nótese que incluso el rmse encontrado es más bajo que con el caso anterior. Por último realizamos el entrenamiento final:

best_hiper <- select_best(grid_fit, metric = "rmse")

final_model_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = best_hiper
                    ) %>%
                    fit(
                      data = data_train
                    ) %>%
                    pull_workflow_fit()

Y evaluamos las predicciones

predicciones <- final_model_tree_fit  %>%
                predict(
                  new_data = data_test_prep,
                  #new_data = bake(transformer_fit, datos_test),
                  type = "numeric"
                )
predicciones %>% head()

Y por último encontramos el error de test:

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

rmse(predicciones, truth = price, estimate = .pred, na_rm = TRUE)

Podemos llegar al mismo reslultado anterior de 61536 con un coste computacional menor.

Modelos

En esta sección veremos modelos disponibles en la paquetería parsnip y al final comparamos sus resultados, estos modelos son:

• GLM (Generalized Linear Models),

• Random forest,

• SVM (Suport Vector Machine), y

• MARS (Multivariate adaptive regression splines).

GLM

Realizaremos el mismo procedimiento que hicimos con árboles de decisión, pero ahora cambiando al modelo de regresión lineal:

# Definición del modelo y los hiperparámetros (GLM)
model_glm <- linear_reg(
                 mode    = "regression",
                 penalty = tune(),
                 mixture = tune()
              ) %>%
              set_engine(engine = "glmnet", nlambda = 100)

# Procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Validación y generación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(model_glm)

# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  penalty(range = c(0, 1), trans = NULL),
                  mixture(range = c(0, 1), trans = NULL),
                  # Número de combinaciones totales
                  levels = c(10, 10)
                )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

Mostramos los mejores resultados de la optimización:

show_best(grid_fit, metric = "rmse", n = 10)

Usamos los hiperparámetros óptimos para obtener el modelo final:

best_hiper <- select_best(grid_fit, metric = "rmse")

modelo_glm <- finalize_workflow(
                x = workflow_modelado,
                parameters = best_hiper
              )

modelo_glm_fit <-  modelo_glm %>%
                   fit(
                     data = data_train
                   )

Calculamos las predicciones con el test

predicciones <- modelo_glm_fit %>%
                predict(
                  new_data = data_test,
                  type = "numeric"
                )

Y por último calculamos el error

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

error_test_glm  <- rmse(
                     data     = predicciones,
                     truth    = price,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "GLM"
                   )
error_test_glm

Random Forest

Un modelo Random Forest está formado por un conjunto de árboles de decisión individuales, cada uno ajustado empleando una muestra bootstrapping de los datos de entrenamiento.

En el entrenamiento de cada árbol, las observaciones se van distribuyendo por nodos generando la estructura del árbol hasta alcanzar un nodo terminal. Cuando se quiere predecir una nueva observación, esta recorre el árbol acorde al valor de sus predictores hasta alcanzar uno de los nodos terminales. La predicción del árbol es la media de la variable respuesta (la moda en problemas de clasificación) de todas las observaciones de entrenamiento que están en este mismo nodo terminal. La predicción de un modelo Random Forest es la media de las predicciones de todos los árboles que lo forman.

Encontremos el modelo utilizando random forest:

# Definición del modelo y los hiperparámetros a optimizar
modelo_rf <- rand_forest(
                 mode  = "regression",
                 mtry  = tune(),
                 trees = tune(),
                 min_n = tune()
              ) %>%
              set_engine(engine = "ranger")

# Procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_rf)

# Grid de hiperparámetros
hiperpar_grid <- grid_max_entropy(
                  # Rango de búsqueda para cada hiperparámetro
                  mtry(range = c(1L, 10L), trans = NULL),
                  trees(range = c(500L, 3000L), trans = NULL),
                  min_n(range = c(2L, 100L), trans = NULL),
                  # Número de combinaciones totales
                  size = 100
                )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

Mostramos los resultados de optimización:

show_best(grid_fit, metric = "rmse")

Usamos los mejores hiperparámetros en el modelo:

best_hiper <- select_best(grid_fit, metric = "rmse")

modelo_rf <- finalize_workflow(
                x = workflow_modelado,
                parameters = best_hiper
             )

modelo_rf_fit <- modelo_rf %>%
                 fit(
                  data = data_train
                 )

Evaluamos las predicciones de test

predicciones <- modelo_rf_fit %>%
                predict(
                  new_data = data_test,
                  type     = "numeric"
                )

y encontramos el error

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

error_test_rf  <- rmse(
                     data     = predicciones,
                     truth    = price,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "RF"
                   )
error_test_rf

SVM

El método de clasificación-regresión Suport Vector Machine (SVM) es un método de no paramétrico, se fundamenta en el concepto de hiperplano y de máxima separación. Los modelos SVM encuentran el hiperplano óptimo capaz de predecir el valor/clasificación de las observaciones con el menor error posible.

Para que el modelo nos funcione correctamente debemos instalar el paquete kernlab.

# Definición del modelo y los hiperparámetros a optimizar
modelo_svm <- svm_rbf(
                 mode      = "regression",
                 cost      = tune(),
                 rbf_sigma = tune(),
                 margin    = tune()
              ) %>%
              set_engine(engine = "kernlab")

# Procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_svm)

# Grid hiperparámetros
hiperpar_grid <- grid_random(
                  # Rango de búsqueda para cada hiperparámetro
                  cost(range = c(-10, -1), trans = log2_trans()),
                  rbf_sigma(range = c(-10, 0), trans = log10_trans()),
                  svm_margin(range = c(0, 0.2), trans = NULL), 
                  # Número de combinaciones totales
                  size = 100
                )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

Mostramos los resultados de la optimización

show_best(grid_fit, metric = "rmse", n = 10)

Y utilizamos los hiperparámetros óptimos para el modelo

best_hiper <- select_best(grid_fit, metric = "rmse")

modelo_svm <- finalize_workflow(
                  x = workflow_modelado,
                  parameters = best_hiper
              )

modelo_svm_fit <- modelo_svm %>%
                  fit(
                    data = data_train
                  )

Encontramos las predicciones en el test:

predicciones <- modelo_svm_fit %>%
                predict(
                  new_data = data_test,
                  type     = "numeric"
                )

Y finalmente encontramos el error

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

error_test_svm <- rmse(
                     data     = predicciones,
                     truth    = price,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "SVM"
                   )
error_test_svm

MARS

Multivariate adaptive regression splines (MARS) es una extensión no paramétrica de los modelos de regresión lineal que automatiza la incorporación de relaciones no lineales entre los predictores y la variable respuesta, así como interacciones entre los predictores.

Para que este funcione bien debemos instalar el paquete earth.

# Definición del modelo y los hiperparámetros a optimizar
modelo_mars <- mars(
                 mode  = "regression",
                 num_terms = tune(),
                 prod_degree = 2,
                 prune_method = "none"
               ) %>%
               set_engine(engine = "earth")

# Procesado
transformer <- recipe(
                  formula = price ~ .,
                  data =  data_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              strata  = price
             )

# Workflow
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_mars)

# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  num_terms(range = c(1, 20), trans = NULL),
                  levels = 20
                )

# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

Mostramos los resultados de la optimización

show_best(grid_fit, metric = "rmse", n = 10)

Utilizamos en el modelo los hiperparámetros óptimos:

best_hiper <- select_best(grid_fit, metric = "rmse")

modelo_mars <- finalize_workflow(
                  x = workflow_modelado,
                  parameters = best_hiper
               )

mmodelo_mars_fit <- modelo_svm %>%
                    fit(
                      data = data_train
                    )

Encontramos las predicciones test:

predicciones <- mmodelo_mars_fit %>%
                predict(
                  new_data = data_test,
                  type     = "numeric"
                )

Encontramos el error

predicciones <- predicciones %>% 
                bind_cols(data_test_prep %>% select(price))

error_test_mars <- rmse(
                      data     = predicciones,
                      truth    = price,
                      estimate = .pred,
                      na_rm    = TRUE
                   ) %>%
                    mutate(
                      modelo = "MARS"
                    )
error_test_mars

Comparación

Error de validación cruzada

set.seed(1234)
cv_folds <- vfold_cv(
              data    = data_train,
              v       = 5,
              repeats = 5,
              strata  = price
            )


registerDoParallel(cores = parallel::detectCores() - 1)

validacion_glm <- fit_resamples(
                    object       = modelo_glm,
                    # preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse),
                    control      = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "GLM")

validacion_svm <- fit_resamples(
                    object        = modelo_svm,
                    #preprocessor = transformer,
                    resamples     = cv_folds,
                    metrics       = metric_set(rmse),
                    control       = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "SVM")

validacion_rf <- fit_resamples(
                    object       = modelo_rf,
                    # preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse),
                    control      = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "RF")

validacion_mars <- fit_resamples(
                     object       = modelo_mars,
                     # preprocessor = transformer,
                     resamples    = cv_folds,
                     metrics      = metric_set(rmse),
                     control      = control_resamples(save_pred = FALSE)
                   ) %>%
                   collect_metrics(summarize = FALSE) %>%
                   mutate(modelo = "MARS")

stopImplicitCluster()

bind_rows(
  validacion_glm,
  validacion_rf,
  validacion_svm,
  validacion_mars
) %>%
ggplot(aes(x = modelo, y = .estimate, color = modelo)) +
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.2) +
  geom_point(alpHa = 0.1) +
  labs(title = "Error de validación cruzada", y = "RMSE") +
  theme_bw() +
  theme(legend.position = "none")

Warning: Ignoring unknown parameters: alpHa

Error de test

errores_test <- bind_rows(
                  error_test_glm,
                  error_test_rf,
                  error_test_svm,
                  error_test_mars
                )

errores_test %>% select(modelo, .metric, .estimate)

ggplot(data = errores_test) +
geom_col(aes(x = modelo, y = .estimate, fill= modelo), color = "gray") +
coord_flip() +
labs(title = "Error de test") +
theme_bw()

Referencias

• Kuhn et al., (2020). Tidymodels: a collection of packages for modeling and machine learning using tidyverse principles. https://www.tidymodels.org

• G. James, D. Witten, T. Hastie, R. Tibshirani. An Introduction to Statistical Learning. MIT Press, 2010.

• Linear Models with R By Julian J. Faraway

• Machine Learning con R y tidymodels, Joaquín Amat Rodrigo.